function [M,V]=initmv(A,B)
% [M,V]=INITMV(A,B)
% calculates initial mean and variance for the shock processes
%
% MODEL
%	x(t+1) = A*x(t) + B*e(t+1)
%	e(t)   ~ iid N(0,I)
% INPUT
%	c, A, B as above
% OUTPUT
%	M: initial/unconditional mean 
%   V: initial/unconditional variance

% Sungbae An
% Updated: 08/17/2005 

n = size(A,1);							% dimension of state
M = zeros(n,1);

% try dlyap.m in Control Toolbox
%if exist('dlyap.m')==2;
try
	% discrete Lyapunov equation solver
	V = dlyap(A,B*B');
catch
	% use Riccatti equation solver
	dum1 = kron(A,A);
	dum2 = B*B';

	V  = reshape(inv(eye(n^2)-dum1)*vec(dum2),n,n);
	V  = symmat(V);								% initial variance

	% test accuracy of the solution of Riccatti equation 
	% if the solution is not accurate enough, solve it numerically.
	tst = riccati(ltrvec(V),dum1,dum2);
	tst = sqrt(tst'*tst);
	if tst > 1e-8
		Vn = fsolve(@riccati,ltrvec(V),optimset('Display','off'),dum1,dum2);
		V  = ltrmat(Vn);
	end
end


function y = riccati(x,a,b)
% RICCATI(X,A,B)
%   evaluates the residual from Riccati equation

% 02/28/2002
% Sungbae An

z = ltrmat(x);
n = size(z,1);   
y = ltrvec(reshape( vec(z) - a*vec(z) - vec(b) ,n,n) );
